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SUMMARY 

We outline specific features of numerical simulations of metamaterial wedges and interfaces. We 
study the effect of different positioning of a grid in the Yee method, which is necessary to obtain 
consistent convergence in modeling of interfaces with metamaterials characterized by negative 
dielectric permittivity and negative magnetic permeability. We demonstrate however that, in the 
framework of the continuous-medium approximation, wave scattering on the wedge may result in a 
resonant excitation of surface waves with infinitely large spatial frequencies, leading to non-convergence 
of the simulation results that depend on the discretization step. Copyright © 2000 John Wiley & 
Sons, Ltd. 
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1. INTRODUCTION 

Recent theoretical 121 El E] and experimental [3 El d studies have shown the possibihty of 
creating novel types of microstructured materials that demonstrate the property of negative 
refraction. In particular, the composite materials created by arrays of wires and split-ring 
resonators were shown to possess both negative real parts of magnetic permeability and 
dielectric permittivity for microwaves. These materials are often referred to as left-handed 
metamaterials, double-negative materials, negative-index materials, or materials with negative 
refraction. Properties of such materials were first analyzed theoretically by V. Veselago a 
long time ago [HI, but they were demonstrated experimentally only recently. As was shown 
by Veselago left-handed metamaterials possess a number of peculiar properties, including 
negative refraction for interface scattering, inverse light pressure, reverse Doppler effect, etc. 

Many suggested and demonstrated applications of negative-index metamaterials utilize 
unusual features of wave propagation and scattering at the interfaces. In particular, the effect 
of negative refraction can be used to realize focusing with a flat slab, the so-called planar 
lens [HI; in a sharp contrast with the well-known properties of conventional lenses with a 
positive refractive index where curved surfaces are needed to form an image. Moreover, the 
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resolution of the negative-index flat lens can be better than a wavelength due to the effect of 
amplification of evanescent modes 

Direct numerical simulations provide the unique means for a design of microwave and 
optical devices based on the negative-index materials, however any realistic simulation should 
take into account metamaterial dispersion and losses |1(J[ 1111 112L [T^ as well as a nonlinear 
response Such numerical simulations are often carried out within the framework of the 
effective medium approximation, when the metamaterial is characterized by the effective 
dielectric permittivity and magnetic permeability. This simplification allows for modelling 
of large-scale wave dynamics using the well-known finite-difference time-domain (FDTD) 
numerical methods jl5| . 

In this paper, we discuss the main features and major difficulties in applying the standard 
FDTD numerical schemes for simulating wave scattering by wedges and interfaces of finite- 
extend negative-index metamaterials, including a key issue of positioning of a discretization 
grid in the numerical Yee scheme |lf)) necessary to obtain consistent convergence in modeling 
surface waves at an interface between conventional dielectric and metamaterial with negative 
dielectric permittivity and negative magnetic permeability. In particular, we demonstrate that, 
in the framework of the continuous-medium approximation, wave scattering on the wedge may 
result in a resonant excitation of surface waves with infinitely large spatial frequencies, leading 
to non-convergence of the simulation results that depend on the discretization step. 



2. BASIC EQUATIONS 

We consider a two-dimensional problem for the propagation of TE-polarizcd electromagnetic 
waves in the plane {x,z), where the medium properties are isotropic and characterized 
by the dielectric permittivity e and magnetic permeability /x. In the absence of losses, 
2m e = Im /i = 0. The response of negative-index materials is known to be strongly frequency 
dependent , however in the linear regime the wave propagation at different wavelengths can 
be described independently. The stationary form of Maxwell's equation for the complex wave 
envelopes is well-known 



dHz dHx iuj dEy iuj dEy iuj 

7^ = —£(x,z)Ey, = fi{x,z)H^, — — = — / 

ox oz c oz c ox c 



where H^, H^, and Ey are the components of the magnetic and electric fields, respectively, lo is 
angular frequency, and c is the speed of light in vacuum. The system of coupled equations (^) 
can be reduced to a single Helmholtz-type equation for the electric field envelope. 



d / 1 dE„\ , B f 1 dE„\ 



where •n?{x, z) = e(x, z)^{x, z), and n is the refractive index of the medium. 



3. WAVE SCATTERING BY A NEGATIVE-INDEX SLAB 

The concept of perfect sub- wavelength imaging of a point source through rcconstitution of the 
evanescent waves by a fiat lens has remained highly controversial (17| because it is severely 
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Figure 1. (a) Geometry of the scattering problem, the slab of the thickness d2 is made of a negative- 
index metamaterial with both 62 and fi2 negative, (b) Discretization scheme in the Yee method. 



constrained by anisotropy and losses of the metamaterials. Nevertheless, several numerical 
studies showed that nearly-perfect imaging should be expected even under realistic conditions 
when both dispersion and losses are taken into account 10, 11, 12, 13 . In this section, we 
consider the numerical simulations of the wave scattering by a slab of the negative-index 
material, i.e. the problem close to that of the perfect lens imaging, and discuss the convergence 
of the Yee numerical discretization scheme. 



3.1. Geometry and discretization 

We start our analysis by considering wave propagation through a slab made of the negative- 
index material, as schematically illustrated in Fig. ^ with homogeneous properties in the 
{y,z) plane characterized by two functions, the electric permittivity s = e{x) and magnetic 
permeability = ^i{x). To solve this problem numerically, we employ the well-known numerical 
Yee method jl5| and perform the discretization of the electric and magnetic fields on a square 
grid {xm = h Zn ~ h n) presenting the fields in the form, 

Ey\m,n — {Ey)m,n, H z\m+l /2,n ^ {H z)\m+l /2,n, Hx\m,n+l/2 = {Hx)\m,n+l/2i (3) 

where we use the notation 

(•)n,m = / / •dxdz. 

Then, we replace the continuum model by a closed set of the discrete equations for the field 
amplitudes obtained by averaging equations ^ over the cells of discretization mesh, taking 
into account the continuity of the tangential field components at the interface |15| . 



Hz\m+l/2,n ~ Hz\m-l/2,n Hx\m,n+l/2 ^ Hx\m.n-l/2 _ i^ 



h c 



Ey\m,n+1 ^-i/lm^n, i> _ 17 I ( A\ 
\M /m,n+l/2 — —-tlx\m,n+l/2^ 

Ey\m,n -^iy|m+l.n i^ j \ tt \ 

~ = ~ Wm+l/2,n"-z\m.+ l/2,n- 
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Figure 2. (a) Spectrum of wavenumbers for a negative-index layer: exact (circles) and discrete 
(crosses, for A'^ = 512) solutions, (b) Absolute differences between the exact and discrete values of 
vs. the number of points (A'^) along x for the marked points 1,2,3 in (a). Bottom: Mode profiles marked 
1,2,3 in (a). The computational domain is < a; < 6.4, 6,2 — 1.4 is the width of the negative-index 
layer with £2 = —1.2 and fi2 — —1.5, and ei — fii — 1. The wavenumber in vacuum is normalized as 

LU/C^ 1. 



Whereas the general form of the discrete equations is well known jl5j , we point out a number 
of specific features arising in numerical simulations of the waves scattering at the interfaces 
with the negative-index media. Since the real parts of both e and /i change sign at these 
interfaces, the corresponding averaged values may become small or even vanish for a certain 
layer position with respect to the numerical grid. In this case, Eqs. Q may become (almost) 
singular, leading to poor convergence. In this paper, we suggest that consistent convergence 
can be achieved by artificially shifting the layer boundary with respect to the grid in order 
to ensure that the averaged values do not vanish. This shift will not exceed h/2, assuring 
convergence as the step-size is decreased. 

Because the tangential component Ey of the electric field should be continuous at the 
interface, is seems that a natural choice is to align the boundary position with the grid 
points Xm, where i?y|m,„ is defined, and we use this configuration in the numerical simulations 
presented below. However, we note that such a selection leads to singularities for averaged 
values if Si = —£2 or /zi = — /X2, which coincides with the flat-lens condition. Therefore, it 
is necessary to take into account losses in the metamaterial, described by nonzero imaginary 
parts of the complex values £2 and fj,2 , or to choose a different boundary alignment to the grid 
for the numerical simulations of perfect lenses jj. 

3.2. Wave spectrum and convergence of discrete solutions 

In order to illustrate the convergence of the proposed numerical scheme, we compare the 
solutions of discrete and continuous equations. We note that wave scattering from an infinite 
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Figure 3. (a) Geometry of the scattering problem, the finite-extent slab is made of a negative-index 
metamaterials with both £2 and jj,2 negative, (b) Discretization scheme in the Yee method. 
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Figure 4. Dependence of (a) the cosine of the singularity parameter and (b) its real (solid) and 
imaginary (dashed) parts on /X2/M1 according to Eq. ®. Shading marks the region with Irwy 7^ 0. 



layer is fully characterized by the properties of spatial modes, which wavevector components 
along the layer (kz) are conserved. These modes have the form 

^{x, z) ~ £(x] kz)exp(ikzz), ll(x, z) ~ Ti.{x; kz)exp{ikzz). (5) 

Substituting Eqs. © into Eq. and Eq. Q), we obtain a set of corresponding continuous 
and discrete eigenmode equations. For every fc^, the mode profiles can be determined 
analytically, e.g., using the transfer-matrix method |18| . The wave spectrum can contain 
solutions corresponding to the guided modes of a negative- index layer and extended 
(or propagating) modes that should also be taken into account as well, in order to describe 
scattering of arbitrary fields. 

We solve the discrete eigenmode equations numerically for the slab geometry with periodic 
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Figure 5. (a) Amplitudes of reflected plane waves vs. the number of points; (b,c) Electric field profiles 
for N = 64 and A'^ — 256. Parameters are the same as in Fig.|5| except ^2 = —3.5. 



boundary conditions, and compare the spectrum of eigenvalues with exact solutions of the 
continuous model. In Fig. [2Ia) , we show a part of the spectrum of the discrete eigenvalues 
(crosses) , which indeed coincides with the exact values (circles) . The rate of convergence can 
be judged from Fig.j^Jb), where the differences between the approximate and exact solutions 
are shown in logarithmic scale. 



4. WAVE SCATTERING BY A WEDGE OF NEGATIVE-INDEX MATERIAL 

One of the fundamental problems in the theory of negative-index metamaterials is the wave 
scattering by wedges (2D|, where convergence of numerical methods can be slow due to the 
appearance of singularities |21| . In this section, we demonstrate that the nature of such 
singularities has to be taken into account when performing FDTD numerical simulations. 

Singularity parameter 

The behavior of the electric and magnetic fields at the wedges between homogeneous materials 
characterized by different values of e and fi was described analytically in the pioneering paper 
of Meixner and further refined in the subsequent studies (see, e.g., Ref. and references 
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Figure 6. (a) Amplitudes of reflected plane waves vs. the number of points; (b,c) Electric field profiles 
for A'' = 64 and A'' = 256. Parameters are the same as in Fig. |5| 



therein). In the case of the TE wave scattering by a negative- index wedge, as schematically 
illustrated in Fig.|3| the amplitudes of magnetic fields exhibit singular behavior at the wedge of 
the order of p'''"^, where p is the distance from the wedge. For a 7r/2 wedge angle, corresponding 
to a corner of a rectangular slab, the coefficient 7 is found as [21] 

^^l--^ = ±2cos(77r/2), (6) 

Ail + ^2 

where ni and ^2 are magnetic permeabilities of the two neighboring media. In the case of 
conventional dielectric or magnetic media with /ij > 0, Eq. ^ has solutions with real 7. 
However, when changes its sign at the interface with a negative-index medium, the coefficient 
7 becomes complex playing role of a singularity parameter. This happens when 

-3/11 < p.2 < -Ml/3, (7) 

so that l/ii - M2I/IM1 + M2I > 2, see Fig.gl 

For real 7 (taking the solution with < 7 < 1, according to Ref. HU), the field amplitudes 
decay monotonously away from the corner. In this case, numerical simulations may be based 
on the simplest discretization, although the convergence rate can be improved by taking into 
account the singular behavior in the discrete equations. 
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Figure 7. Same as in Fig.|Slbut for /i = —1.5 + i 10 ^ . 

However, the case of complex 7 corresponds to the fields that oscillate infinitely fast near 
the corner because 

^7-1 ^ pTCe(7-l) log(p)] . 

The second multiplier indicates the excitation of infinitely large spatial harmonics, however 
such a situation is unphysical because the effective-medium approximation of the negative- 
index materials is valid for slowly varying fields only. Therefore, in numerical simulations it is 
necessary to take into account the physical effects that suppress such oscillations, in particular 
we discuss the effect of losses in Sec. 14.31 below. 

4-2. Numerical results 

We now analyze convergence of the numerical finite-difference solutions for the problem of wave 
scattering by a finite-extent negative-index slab. We align the boundaries of the negative-index 
domain with the grid, as shown schematically in Fig.QJb). Since the electric field components 
Ey are continuous at the interfaces, it is possible to obtain the discrete equations that have 
the form of Eqs. 0). 

In order to construct the full solutions for scattering problem, we decompose the field into a 
set of eigenmodes of the negative-index layer (z > 0) and free space (z < 0). More specifically. 



Copyright © 2000 John Wiley & Sons, Ltd. 
Prepared using jnmauth.cis 



Int. J. Numer. Model. 2000; 00:1-6 



WAVE SCATTERING BY METAMATERIAL WEDGES AND INTERFACES 



9 



for z > we have 

EU.n = ^ A,5U^„(fc(^"))exp(ifc(^")z„), HU,„ = ^A,WU^„(fcp"))exp(ifcp')z„), (8) 

j j 

where j is the number of eigcnmodes. Here the summation is performed over the propagating 
modes [Imkz = 0) which transport energy away from the interface, and evanescently decaying 
modes with Xmkz < 0. In free space at z < 0, the field is composed of incident and reflected 
plane waves, 

N/2 N/2 

Ey\rn,n^ ^ Fj''^hxp{i27Tjm/N + iK^^hn) + ^ f]'""''^ exp(i27rjm/7V - ix(-'')z„). 

j=-N/2+l j=~N/2+l 

(9) 

Here N is the number of points in the x direction, and the discrete wavenumber is K'^^^ = 

1 /2 

±sin-i [h^iuyAc^ - sin^jn/N)] ' 2/h The sign of K^^^ is chosen with a proper wave 
asymptotic behavior, i.e., we choose ne{K'-jy) > 0, if K'-J^ is real, and Im{K^^^) > if K<^^^ 
is complex. The magnetic field at z < is found from Eq. Q), with homogeneous parameters 
(e) = £i, (fi) = ^1, (^^^) = ^ii^ ■ Then, we substitute Eqs. ©, © into the first of Eq. |0J, 
and using the condition of the continuity of the electric field, obtain a set of equations for all 
m — 1, . . . ,N that are used to calculate the amplitudes pj^'^^^^ and Aj of the reflected and 
transmitted waves, 

N/2 

j 3 = -N/2+l 

^x|™,-l/2(fcP^)+H,U,i/2(4'^) 



(10) 



9 

(_pM „^j>^<=fl))sin(ifO')/i)exp(i27rjm/iV). 



These equations are solved using the standard linear algebra package. 

We now present results of our numerical simulations for the scattering of normally incident 
plane waves, with Fg™'' = 1, Fj^^ = 0. First, we consider scattering by a negative-index slab 
with /i2 < —3. In this case, we observe a steady convergence of numerical solutions, as shown 
in Fig. [SI This demonstrates that even simplest finite-difference numerical schemes can be 
successfully employed to model the scattering process when the sinularity parameter 7 is real, 
is in a full agreement with earlier studies of wave scattering at dielectric wedges ^T? . 

However, the situation changes dramatically when 7 is complex, i.e. for ^2 — —1.5. According 
to the analytical solution, in this case the magnetic field should oscillate infinitely fast in 
the vicinity of the corner, corresponding to excitation of infinitely large spatial harmonics. 
However, such behavior cannot be described by discrete equations, and we find that in this 
regime solutions of finite- difference equation do not converge, as demonstrated in Fig.|H| 

4-3. Effects of losses 

The analytical description of the edge singularities discussed above is only valid for lossless 
media, i.e. when all e and fi are real. However, the negative-index metamaterials always have 
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Figure 8. Same as in Fig.Hbut for /i = -1.01 + i 10" 



non-vanishing losses, and we have studied whether this important physical effect can regularize 
the field oscillations at the corner. However, our results demonstrate that even substantial 
losses may not be sufficient enough to suppress such oscillations, as presented in the example 
of Fig.d 



4-. 4- Singularities and perfect lenses 

Finally, we consider the problem of wave scattering from the corners of perfect lenses, where 
TZe (£2) — —1 and TZe (^2) — —1 (we take ei = /xi = 1). This is a special case, where the type 
of singularity becomes indefinite if losses are neglected. We find that introducing sufficiently 
large losses does indeed regularize the field oscillations at the corners , leading to convergence 
of numerical simulations, as demonstrated for the example of Fig.|Hl However, this only occurs 
when the value of losses exceeds a certain threshold; if the losses are too weak then non- 
convergent behavior is again observed, as shown in Fig. El The threshold value of losses for 
achieving convergence of the numerical scheme is increased for larger \TZe {^2) + 1|- We note 
that this is completely different from the temporal dynamics at an infinitely extended slab, 
where convergence to a steady state is eventually achieved with arbitrarily small losses jl2| . 
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Figure 9. Same as in Fig. Hbut for /i = -1.01 + i 10"^. 
5. CONCLUSIONS 

We have discussed a number of specific features manifested in numerical simulations of wedges 
and interfaces of metamaterials, i.e. composite materials with negative dielectric permittivity 
and negative magnetic permeability. We have demonstrated that a numerical discretization 
grid in the Yee method may have a dramatic effect on the convergence in numerical modelling 
of surface waves at interfaces and wedges. In the framework of the continuous-medium 
approximation, wave scattering on the wedge may result in a resonant excitation of surface 
waves with infinitely large spatial frequencies, leading to non-convergence of the numerical 
simulation results that depend strongly on the value of the discretization step. We find that 
sufficiently high losses may suppress oscillations and allow to obtain converging solutions to 
the scattering problem, however in the case of smaller losses it may be necessary to take into 
account the meta-material properties beyond the effective-medium approximation, such as the 
effect of spatial dispersion. 
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